Анализ временного ряда

Временной ряд

Данные – траты жителей Великобритании за рубежом с января 1980г. по декабрь 2016г. в миллионах £.

Чуть более свежие данные можно найти здесь: https://www.ons.gov.uk/peoplepopulationandcommunity/leisureandtourism/timeseries/gmam/ott

uk <- read.csv("UK.csv", sep = ',', stringsAsFactors = FALSE)
uk <- ts(uk$EXPEND, start = c(1980, 1), frequency = 12)
df.uk <- data.frame(seq.Date(from = as.Date("1980-01-01"), to = as.Date("2016-12-01"), by = "month"), uk)
colnames(df.uk) <- c("DATE", "EXP")
tail(uk, 24)

      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
2015 2371 2071 2548 2901 3255 3690 4442 5506 4461 3445 2365 1972
2016 2448 2523 2607 3529 3489 3969 4569 6097 5215 4030 2820 2120

plot(uk)

Тренд сначала похож на линейный, примерно в 2008 году резко меняется. Имеет смысл рассматривать только конец ряда, сославшись на то, что резкое изменение тренда связано с финансовым кризисом.

Модель не аддитивная, амплитуда сезонной компоненты не постоянна. Посмотрим, что происходит при логарифмировании.

plot(log(uk))

Амплитуда не постоянна, значит, модель не является мультипликативной.

Периодограмма

\((x_1, \ldots, x_N)\) – ряд.

Представим его в виде \(x_n = C_0 + \sum\limits_{k=1}^{[(N-1)/2]}(C_kcos(2\pi n k/N) + S_ksin(2\pi n k/N))\) (\(+ C_{N/2}(-1)^n\)), если \(N\) – четное.

Общий вид периодограммы:

\(\widetilde{П}_x^N(w) = \frac{1}{N}|\sum\limits_{n = 0}^{N - 1} e^{-2\pi iwn} x_{n + 1}|^2\), \(w \in (-1/2, 1/2)\)

Так как ряд вещественный, рассматриваем \(П_x^N(k/N) = \frac{N}{2}\begin{cases} 2C_0^2, &k = 0\\ C_k^2 + S_k^2, &0 < k < \frac{N}{2}\\ 2C_{N/2}^2, &k = \frac{N}{2} \end{cases}\)

spec.pgram(uk, detrend = FALSE, fast = FALSE, log='no', xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

spec.pgram(uk, detrend = TRUE, fast = FALSE, log='no', xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

Заметны периоды 12 (год), 6 (полгода), 4 (квартал), 2.4, 2 и частота близкая к нулю.

Отступление: Пример про растекание частоты

\(x_n = cos(2\pi n/10)\)

par(mfrow=c(1,2))
spec.pgram(cos(2*pi*1:92/10), detrend = FALSE, fast = FALSE, log='yes', taper=0, pad=0) #растекается
spec.pgram(cos(2*pi*1:100/10), detrend = FALSE, fast = FALSE, log='yes', taper=0, pad = 0)#не растекается

Отступление: белый шум

• \(w_t, t = 0, \pm 1, \pm 2, \ldots\) – некоррелированые случайные величины, \(\mathbb{E}w_t = 0,~\sigma_w^2\), пишут, например, \(w_t \sim wn(0,~\sigma_w^2)\).

• если кроме того \(w_t\) – iid, пишут \(iid(0,~\sigma_w^2)\)

• если \(w_t\) – iid и \(w_t \sim N(0, ~\sigma_w^2)\), шум гауссовский

Периодограмма гауссовского шума

wn <- ts(rnorm(1000))
spec.pgram(wn, detrend = FALSE, log = 'no')

spec.pgram(wn, detrend = FALSE, log = 'yes')

Автоковариационная функция

acf(wn)

Cглаживание временного ряда и выделение тренда

Тренд – медленно меняющаяся компонента ряда.

Тренд – случайный процесс, где медленно меняющаяся компонента превалирует.

Чем отличается сглаживание от выделения тренда?

• Сглаживание даст медленно меняющуюся компоненту. Если понимать тренд как некоторую неслучайную функцию, то можем получить тренд, а можем и нет.

Линейный фильтр, FIR, причинный фильтр, АЧХ, ФЧХ

\(x_n\) – ряд.

Фильтры

• \(y_j = (\Phi(X))_j = \sum\limits_{i = -\infty}^{\infty} h_i x_{j- i}\) – линейный фильтр

• \(\{h_i\}\) – импульсная характеристика (impulse response)

• FIR – finite impulse response, конечное число \(h_i\) ненулевые, т.е. \((\Phi(X))_j = \sum\limits_{i = -r_1}^{r_2} h_i x_{j- i}\)

• Причинный фильтр (casual filter), смотрим только в “прошлое”: \((\Phi(X))_j = \sum\limits_{i = 0}^{r} h_i x_{j- i}\)

АЧХ, ФЧХ

• \(H_{\Phi}(z) = \sum\limits_i h_i z^{-i}\) – передаточная функция (transfer function)

• \(|H_{\Phi}(e^{i2\pi w})|\) – Амплитудно частотная характеристика (АЧХ, frequency response), показывает как меняется амплитуда в зависимости от частоты

• \(\varphi_\Phi(w) = Arg H_{\Phi}(\varphi^{i2\pi w})\) – Фазочастотная характеристика (phase response)

Пример

Функция скользящего среднего

MovAvg <- function(series, M){
  avg <- sapply(1:length(series), FUN = function(i){
    if(M %% 2 == 0){
      l <- ((i - M/2 + 1):(i-1))
      r <- (i:(i + M/2))
    } else{
       l <- ((i - floor(M/2)):(i-1))
      r <- ((i + 1):(i + floor(M/2)))
    }
    
    l <- l[l > 0]
    r <- r[r <= length(series)] 
    h <- 1/(length(r) + length(l) + 1*(M %% 2))
    if(M %% 2 == 0) sum(series[c(l,r)])*h 
    else sum(series[c(l, i, r)])*h
  })
}

Скользящее среднее по прошлому

MovAvg.cause <- function(series, M){
  avg <- sapply(1:length(series), FUN = function(i){
    l <- ((i - M + 1):i)
    l <- l[l > 0]
    h <- 1/length(l)
    sum(series[l])*h 
  })
}

Функция АЧХ

frequencyResponse <- function(h, w){
  fr <- w
  if(max(fr) > 0.5) fr <- fr / (max(w) * 2)
  sapply(fr, function(freq){
    cs <- exp(-(1:length(h))*1i * 2*pi* freq)
    pr <- h * cs
    return(abs(sum(pr)))
  })
}

Пусть \(x_n = cos(2\pi\omega n)\).

Применим фильтр, получим \((\Phi(X))_n = A_{\phi}(\omega)cos(2\pi\omega n + \varphi_{\phi}(\omega))\). АЧХ фильтра – то как меняется амплитуда у косинуса с частотой \(\omega\).

В случае \(x_n = cos(2\pi\omega n)\), \(A_{\phi}(\omega) = 2sin(\pi\omega)\)

xn <- cos(2*pi* 1/12 * (1:100) + 0.5)
plot(xn, type = 'l')

Перейдем к разностям

xn.filtered <- diff(xn)
plot(xn.filtered, type = 'l')

АЧХ

h <-  c(-1,1)
o <-  seq(0, pi, 0.001)/2/pi
plot(frequencyResponse(h,o) ~ o, type = 'l', xlab = "frequency", ylab = "response", main = "Frequency Response")

\(A_{\phi}(\omega) = 2sin(\pi\omega)\)

plot(2*sin(pi*o) ~ o, type = 'l', xlab = "frequency", ylab = "response", main = "Frequency Response")

Роль нормы коэффициентов фильтра

• \(x_n\) – некоторый ряд
• Применили фильтр \(y = \sum\limits_{i = 1}^n a_i x_i\)
• Добавили ошибки \(x_i = x_i + \epsilon_i\), \(D(\epsilon_i) = \sigma^2\), \(\epsilon_i\) – независимы
• Тогда \(\widehat{y} = \sum\limits_{i = 1}^n a_i (x_i + \epsilon_i)\)
• \(\mathbf{E}\widehat{y} = y\), \(D\widehat{y} = \sum\limits_{i = 1}^n a_i^2 \sigma^2 = \|A\|^2\sigma^2\)
• Минимизируя норму, минимизируем дисперсию. Лучший результат при \(a_i = 1/n\)

Cкользящее среднее

Скользящее среднее – линейный фильтр, для которого \(h_i = \frac{1}{2M + 1}\).

Функция, которая строит три периодограммы сразу (для исходного ряда, результата применения фильтра и остатка)

periodogramBuilder <- function(series, filtered){
   par(mfrow=c(1,3))
   
   sp.initial <- spec.pgram(series, detrend = FALSE, fast = FALSE, log='no', taper=0)
   sp.filtered <- spec.pgram(na.omit(filtered), detrend = FALSE, fast = FALSE, log='no', taper=0)
   sp.difference <- spec.pgram(na.omit(series - filtered), detrend = FALSE, fast = FALSE, log='no', taper=0)
   par(mfrow=c(1,1))
}

Функция, которая строит все нужные картинки сразу.

movingAverageFiltering <- function(timeSeries,m){
  if( !is.ts(timeSeries)){
    stop("Object given is not a ts object!")
  }
  
  par(mfrow=c(1,1))
  #Initial
  plot(timeSeries, type='l', lty=2, ylab = 'Rate')
  
  h <- rep(1,(m + 1*(m%%2)))/(m + 1*(m%%2))
  
  #average
  filtered <- ts(MovAvg(timeSeries, m), start = c(1980, 1), frequency = 12) 
  lines(filtered, col='red')
  #difference (detrended)
  difference <- timeSeries - filtered
  plot(difference, type = 'l')
  
  #comparing periodograms
  periodogramBuilder(timeSeries, filtered)
  #Frequency response
  par(mfrow=c(1,1))
  sp.filtered <- spec.pgram(na.omit(filtered), detrend = FALSE, fast = FALSE, log='no', taper=0, plot = FALSE)
  fr.filtered <- frequencyResponse(h, sp.filtered$freq)
  plot(x = sp.filtered$freq, y = fr.filtered, type = 'l', xlab = "frequency", ylab = "response", xaxt = 'n')
  axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))
}

Сгладим ряд. Выберем окно кратное 12, чтобы убрать сезонность.

movingAverageFiltering(uk, 12)

Вся периодика ушла в остаток. Компонента, соответствующая низкой частоте, смешалась с трендом.

По АЧХ видно, что частоты, соответствующие сезонности (1/12, 2/12 и т.д.) фильтр обращает в 0.

Увиличиваем окно

movingAverageFiltering(uk, 24)

movingAverageFiltering(uk, 36)

movingAverageFiltering(uk, 60)

Чем больше окно, тем сильнее подавляем высокие частоты, что видно по АЧХ.

Остановимся на окне 12. Каких-то серьезных улучшений бОльшие окна не дают, но зато теряется больше информации на концах ряда (получим запаздываение).

“Low-pass filter” – высокие частоты убираем, низкие оставляем.

Смещение при сглаживании фильтром скользящего среднего. Роль второй производной

Общий вид фильтра: \[\begin{align*} y(a) = \int_{-\delta}^{\delta} f(a + x) w(x) dx, \quad \int_{-\delta}^{\delta} w(x)dx = 1,\quad \int_{-\delta}^{\delta} x w(x)dx = 0. \end{align*}\]

Пусть \(f\) – некоторая гладкая функция, применим к ней фильтр.

\[\begin{align*} f(a + x) = f(a) + f'(a)x + f''(a)\frac{x^2}{2} + \ldots,\\ y(a) \approx f(a) + f'(a) 0 + \frac{f''(a)}{2} \int_{-\delta}^{\delta} x^2 w(x)dx. \end{align*}\]

В случае скользящего среднего \(w(x) = \frac{1}{2\delta}\), следовательно, у нас всегда будет смещение \(\approx \frac{\delta^2}{3} \frac{f''(a)}{2}\).

Отступление: перерисовка и запаздывание

Перерисовка

Возникает из-за того, что нам не хватает точек справа, чтобы посчитать среднее. При добавлении новой точки среднее пересчитывается.

Последовательно добавляем последние 100 точек.

plot(y = uk, x = (1:444), type = 'l')
for(i in 100:1){
  ma <- MovAvg(head(uk, -i), 12)
  lines(ma, type = 'l', col = 'red', ylim = c(0, 6000), xlim = c(0, 450))
}

Запаздывание

Рассмотрим причинный фильтр, например скользящее среднее, которое считается только по “прошлому”.

plot(y = uk, x = (1:444), type = 'l')
ma.c <- MovAvg.cause(uk, 12)
ma <- MovAvg(uk, 12)
lines(ma.c, type = 'l', col = 'red', ylim = c(0, 6000), xlim = c(0, 450))

lines(ma, type = 'l', col = 'blue', ylim = c(0, 6000), xlim = c(0, 450))
legend("topleft", c("Usual", "Casual"), lty=c(1,1), lwd=c(2.5,2.5),col=c("blue","red"))

Видно, что casual фильтр “догоняет” обычное среднее только через некоторое количество точек.

Скользящая медиана

В отличие от скользяещего среднего, устойчива к аутлаерам. Но негладкая и не является оценкой м.о. для несимметричного ряда. Иногда сначала применяют скользящую медиану, а потом результат сглаживают скользящим средним.

movingMedianFiltering <- function(timeSeries,k){
  if( !is.ts(timeSeries)){
    stop("Object given is not a ts object!")
  }
  
  par(mfrow=c(1,1))
  #Initial
  plot(timeSeries, type='l', lty=2, ylab = 'Rate')
  #median
  filtered <- runmed(timeSeries, k)
  lines(ts(filtered, start = c(1980, 1), frequency = 12 ) , col='red', type='l')
  #difference (detrended)
  difference <- timeSeries - filtered
  plot(difference, type = 'l')
  
  #comparing periodograms
  
  periodogramBuilder(timeSeries, filtered)
}
movingMedianFiltering(uk, 13)

Отступление: ряд без сезонности, с аутлаерами

Рассмотрим ряд среднегодовых отклонений температуры

plot(globtemp)

spec.pgram(globtemp, detrend = FALSE, log = 'no')

spec.pgram(globtemp, detrend = TRUE, log = 'no')

Видно, что сезонности нет.

Добавим аутлаеры.

temp.w.outliers <- globtemp
set.seed(1)
out <- base::sample(1:length(globtemp), 10)
temp.w.outliers[out] <- c(3.3,-3,4,-4,5,-5,4,-4, 5, -5)

Применим скользящую медиану c маленьким окном.

  plot(temp.w.outliers)
  temp.med <- runmed(temp.w.outliers, 3)
  lines(ts(temp.med, start = 1880, frequency = 1) , col='red', type='l')

Избавились от аутлаеров.

Разность первого порядка

df <- diff(uk)
h.df <- c(-1,1)
plot(df, type='l', main='difference')
par(mfrow=c(1,2))

sp.initial <- spec.pgram(uk, detrend = FALSE, fast = FALSE, log='no', taper=0, xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

sp.df <- spec.pgram(na.omit(df), detrend = FALSE, fast = FALSE, log='no', taper=0, xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

par(mfrow=c(1,1))

fr.df <- frequencyResponse(h.df, sp.df$freq/12)
plot(x = sp.df$freq, y = fr.df, type = 'l', xlab = 'frequency', ylab = 'response', xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

“High-pass filter” – пропускает высокие частоты, низкие убирает.

Почему переход к разностям это хорошо?

• Если модель нашего ряда имеет вид \(x_t = \mu_t + y_t\), где \(\mu_t\) – тренд, \(y_t\) – стационарный процесс, то при переходе к разностям получим стационарный ряд, особенно, если тренд фиксирован (например, линейный)

• Увеличивем вклад высоких частот. Например, легче заметить сезонность, как на периодограмме выше.

• Уберем линейный тренд (если перейдем к разностям второго порядка, то квадратичный и т.д.)

Минусы

• Увеличиваем вклад шума. Например, на периодограмме выше не ясно относить ли синус с частотой чуть больше 4/12 к сезонности или шуму.

Сглаживание и выделение тренда с помощью регрессии

Линейная регрессия

fit <- lm(df.uk$EXP ~ df.uk$DATE, df.uk)
plot(df.uk, type='l')
abline(fit, col='red')

periodogramBuilder(df.uk$EXP, fit$fitted.values)

plot(fit$residuals, type='l')

Полиномиальная регрессия

cs <- cos(2*pi*(1:length(uk))/12)
sn <- sin(2*pi*(1:length(uk))/12)
fit <- lm(EXP ~ stats::poly(df.uk$DATE, 10)
                  , df.uk)
plot(df.uk, type='l', ylab="EXP")
lines(y= fitted(fit), x= df.uk$DATE, col='red')

plot(fit$residuals, type='l')

periodogramBuilder(df.uk$EXP, fit$fitted.values)

Kernel Smoothing (kernel regression)

Kernel Smoothing – скользящее среднее, в котором используется весовая функция (kernel).

Cчитаем, что наш ряд имеет вид: \[\begin{align*} x_t = f_t + y_t, \end{align*}\] где \(f_t\) – некоторая гладкая функция, \(y_t\) – стационарный процесс. Тогда \[\begin{align*} \widehat{f}_t = \sum_{i = 1}^{n}w_i(t)x_i, \end{align*}\] где \[\begin{align*} w_i(t) = K(\frac{t-i}{b}) / \sum_{j = 1}^{n}K(\frac{t-j}{b}). \end{align*}\]

Обычно \(K(z) = \frac{1}{\sqrt{2\pi}} \exp(-z^2 / 2)\).

plot(uk, type='l', ylab="EXP")
trend.sm <- ksmooth(time(uk), uk, 'normal', bandwidth = 2)
lines(trend.sm, col='red')

plot(uk - trend.sm$y)

periodogramBuilder(uk, trend.sm$y)

Nearest Neighbor Regression

Строим линейную регрессию на \(k\) ближайших соседей. Т.е. предсказываем \(x_t\) по \(\{x_{t-k/2},\ldots, x_t,\ldots, x_{t+k/2}\}\).

plot(uk, type='l', ylab="m £")
neigh.t <- knn.reg(train = time(uk), y = uk, k = 24)$pred#supsmu(time(uk), uk, span = .3)
neigh.t <- ts(neigh.t, start = c(1980, 1), frequency = 12)
lines(neigh.t, col='red')

plot(uk - neigh.t)

periodogramBuilder(uk, neigh.t)

Сравнение методов

• Линейная регрессия – подходит только для выделения линейного тренда
• Полиномиальная регрессия – не подходит для прогнозирования
• Параметрические методы плохо сработают, если не предполагается, что тренд – всегда одна и та же функция
• Nearest Neighbor Regression, Kernel Smoothing – непараметрические методы, нет ограниченности моделью

Для выделения тренда я бы выбрал Nearest Neighbor Regression или Kernel Smoothing.

Seasonal Decomposition

Хотим представить ряд в виде \(X_n = T_n + S_n + N_n\).

Будем делать все для логарифмированного ряда, так как он похож на линейную модель (амплитуда примерно одинакова). Можно заменить далее сложение на умножение, разность на деление и получить алгоритм для мультипликативной модели.

Хорошее описание алгоритма нашел здесь: http://www.abs.gov.au/websitedbs/d3310114.nsf/51c9a3d36edfd0dfca256acb00118404/c890aa8e65957397ca256ce10018c9d8!opendocument

• Начальная оценка тренда. Применяем скользящее среднее с длиной окна T, в нашем случе кратной 12.

T1 <- MovAvg(log(uk), 24)
SN1 <- log(uk) - T1
periodogramBuilder(log(uk), T1)

• Получили ряд в виде \(X_n = \widetilde{T}_n + \widetilde{S_n + N_n}\)

• Оценка сезонной компоненты. Сглаживаем \(\widetilde{S_n + N_n}\) с маленьким окном. Получаем \(\widetilde{S_n + N_n} = \widetilde{S}_n + \widetilde{N}_n\).

• Убираем полученную сезонность из исходного ряда, получаем оценку исправленного ряда (adjusted).

S1 <- MovAvg(SN1, 2)
ADJ1 <- log(uk) - S1
periodogramBuilder(SN1, S1)

• Улучшенная оценка тренда. Применяем скользящее среднее с большим окном к исправленному ряду \(X_n - \widetilde{S}_n\), получаем \(\widetilde{\widetilde{T}}_n + \widetilde{\widetilde{N}}_n\)

• Отсюда второй раз оцениваем \(\widetilde{\widetilde{S_n + N_n}} = X_n - \widetilde{\widetilde{T}}_n\)

T2 <- MovAvg(ADJ1, 24)
SN2 <- log(uk) - T2
periodogramBuilder(ADJ1, T2)

• Применяем скользящее среднее с маленьким окном к \(\widetilde{\widetilde{S_n + N_n}}\), получаем \(\widetilde{\widetilde{S}}_n + \widetilde{\widetilde{\widetilde{N}}}_n\). Cнова оцениваем исправленный ряд

S2 <- MovAvg(SN2, 2)
ADJ2 <- log(uk) - S2
periodogramBuilder(SN2, S2)

• Финальная оценка тренда и шума

T3 <- MovAvg(ADJ2, 24)
N1 <- ADJ2 - T3
periodogramBuilder(ADJ2, T3)

• Получили \(X_n = \widetilde{\widetilde{T}}_n + \widetilde{\widetilde{S}}_n + \widetilde{\widetilde{\widetilde{N}}}_n\)

par(mfrow=c(4,1))
plot(log(uk), type = 'l', ylab='Original')
plot(T3, type='l', ylab='Trend')

plot(S2, type='l', ylab = 'Seasonal Component')
plot(N1, type='l', ylab='Noise')

acf(N1, lag.max = 100)

Шум белым не получился. Как видно из последних периодограмм, много периодических компонент так и не выделились в сезонную компоненту.

STL (LOESS)

Алгоритм

Inner loop

Повторяем \(n_{(i)}\) раз.

• Detrending

\(Y_n - T_n^{(k)}\), \(T_n^{(0)} = 0\), \(Y_n\) – исходный ряд, \(T_n^{(k)}\) – оценка тренда на шаге k внутреннего цикла

• Cycle-subseries smoothing

Каждый период длины \(n_{(p)}\) сглаживается LOESS(q), \(q = n_{(s)}, d = 1\), где \(d\) – степень полинома, \(q\) – число соседей.

Кроме того, считаем сглаженные значения для одной точки перед началом и одной точки после конца каждого периода. Это нужно для того, чтобы после применения скользящих средних на следующем шаге, осталось N точек. То есть, чтобы мы ничего не потеряли.

Получили ряд \(C_{n}^{(k + 1)}\) содержащий \(N + 2*n_{(p)}\) точек.

• Low-Pass Filtering of Smoothed Cycle-Subseries

Последовательно применяем 2 скользящих средних с длинной окна \(n_{(p)}\).

Применяем скользящее среднее с длинной окна 3.

Применяем LOESS с параметрами \(q = n_{(l)}, d = 1\).

Получили ряд \(L_n^{(k + 1)}\)

• Detrending of Smoothed Cycle-Subseries

Считаем ряд \(S_n^{(k + 1)} = C_n^{(k + 1)} - L_n^{(k + 1)}\).

• Deseasonalizing

Считаем ряд \(Y_n - S_n^{(k + 1)}\).

• Trend Smoothing

Ряд без сезонности, вычисленный на предыдущем шаге, сглаживается LOESS с \(q = n_{(t)}, d = 1\).

Получили оценку тренда \(T_n^{(k + 1)}\)

Outer Loop

Повторяем \(n_{(o)}\) раз.

Пусть мы получили оценки \(T_n, S_n\) c помощью внутреннего цикла, тогда остаток \(R_n = Y_n - T_n - S_n\).

Определим веса, которые позволят нам исключить аутлаеры.

\(h = 6 median(|R_n|)\), тогда вес в точке \(n\) \(\rho_n = B(|R_n|/h)\), где \(B\)-биквадратная весовая функция.

На следующей итерации во внутреннем цикле, на шагах 2 и 6, веса, которые присваиваются точкам в LOESS умножим на соответствующие \(\rho_n\).

Параметры

Таким образом, STL имеет 6 параметров:

• \(n_{(p)}\) – количество наблюдений в каждом цикле сезонной компоненты (период)
• \(n_{(i)}\) (inner) – число итераций внутреннего цикла
• \(n_{(o)}\) (outer) – число итераций внешнего цикла
• \(n_{(l)}\) (l.window) – сглаживающий параметр для low-pass фильтра (количество соседей в LOESS)
• \(n_{(t)}\) (t.window) – сглаживающий параметр для тренда (количество соседей в LOESS)
• \(n_{(s)}\) (s.window) – сглаживающий параметр сезонности (количество соседей в LOESS)

Другие парамаетры, доступные в функции

• t(s,l).degree – степень полинома в LOESS для соответствующих компонент (0 или 1)

Использование

Выберем следующие параметры

• s.window = 13 так как у нас годовая периодичность, и параметр должен быть нечетным
• l.window = 13 рекомендуется брать ближайшим нечетным к \(n_{(p)}\)
• outer = 0, аутлаеров нет
• inner = 2, утверждается, что достаточно и одной итерации внутреннего цикла, но для лучшей сходимости рекомендуют 2
• t.window = 23, рекомендуется брать \(n_{(t)} = [1.5 n_{(p)}/(1 - 1.5/n_{(s)})]_{odd}\)

uk.stl <- stl(uk, s.window = 13, l.window = 13, outer = 0, inner = 2, t.window = 23)
plot(uk.stl)

acf(uk.stl$time.series[,3], lag.max = 100)

Стабилизация дисперсии

Остаток после stl

rem <- uk.stl$time.series[,3]
plot(rem)

Прологорифмируем исходный ряд и снова посмотрим на остаток после выделения тренда.

uk.log <- log(uk)
uk.log.stl <- stl(uk.log, s.window = 13, l.window = 13, outer = 0, inner = 2, t.window = 23)
log.rem <- uk.log.stl$time.series[,3]
plot(log.rem)

Оценим дисперсию

Идея * Выделив тренд, получили шум с некоторой дисперсией \(\xi_n = \sigma(n)\epsilon_n\)

• Выделить тренд у \(\xi_n^2 = \sigma(n)^2\epsilon_n^2\) – все равно, что найти среднее \(\mathbf{E}\xi_n^2 = \sigma(n)^2\)

Возведем остаток в квадрат и выделим тренд скользящим средним

dif.sqr <- ts(log.rem^2)
plot(dif.sqr, type = 'l')
sigma.n.sqr <- MovAvg(dif.sqr, 24)
lines(sigma.n.sqr, col='red')

legend("topleft", c("Rem^2", "sigma(n)^2"), lty=c(1,1), lwd=c(2.5,2.5),col=c("black","red"))

Построим огибающую

noise <- ts(log.rem / sqrt(sigma.n.sqr))
plot(log.rem, type = 'l')
lines(ts(sqrt(sigma.n.sqr), start = c(1980, 1), frequency = 12), col = 'blue')

lines(ts(-sqrt(sigma.n.sqr), start = c(1980, 1), frequency = 12), col = 'blue')

Оценка шума, ее периодограмма и автоковариационная функция

plot(noise)

spec.pgram(na.omit(noise), detrend = FALSE, fast = FALSE, log='no', taper=0)

acf(noise, lag.max = 50)